Entropy production in the cyclic lattice Lotka-Volterra model 
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The cycUc Lotka-Volterra model in a D- dimensional regular lattice is considered. Its "nucleus 
growth" mode is analyzed under the scope of Tsallis' entropies Sq — (1 — ^iPi)/{q— 1), ^ G M. It is 
shown both numerically and by means of analytical considerations that a linear increase of entropy 
with time, meaning finite asymptotic entropy rate, is achieved for the entropic index = 1 — 1/^- 
Although the lattice exhibits fractal patterns along its evolution, the characteristic value of q can 
be interpreted in terms of very simple features of the dynamics. 
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I. INTRODUCTION 

Over a decade ago, Tsallis proposed an entropic 
form as a starting point for a possible generalization of 
Boltzmann-Gibbs statistics [1] (see also [2] for a review 
on the subject). The generalized entropy has the form 



So 



k- 



1 



with q G 



(1) 



where /c is a positive constant and p a normalized prob- 
ability density. The usual Boltzmann-Gibbs- Shannon 
(BGS) entropy is recovered when q = \ [Si = 
—k J dxp{x) lnp{x)]. 

Many experimental and numerical data are well ap- 
proximated by g'-exponentials (the probability distribu- 
tions that maximize Sq under simple constraints), giving 
support to the applicability of the new entropic measure 
to those systems. In the case of the high-dimensional 
systems of interest in statistical physics, first-principle 
derivations are still to be made. However, Tsallis' en- 
tropy seems to be the appropriate one for certain low- 
dimensional dynamical systems, such as unimodal maps 
at the edge of chaos, where the usual Lyapunov exponent 
vanishes. In such cases the characteristic value of q is re- 
lated to relevant properties of the dynamics and it can 
be calculated by several paths, such as the multifractal 
spectrum of the attractor [3] or the sensitivity to initial 
conditions [4], or even, from a rigorous renormalization 
group analysis [5]. Furthermore, this scenario has been 
shown to be valid also for two-dimensional maps [6]. 

One of the properties of Sq that, for an appropriate 
value of the entropic index, can make it in some cases 
preferable to the standard entropy is the possibility of 
having an asymptotic rate of entropy production with a 
non-trivial value. Among the a priori infinite possible 
values of g, the one leading to a linear increase of the en- 
tropy with time, implying finite rate of entropy growth, is 
selected as being qc^ a value of q characteristic of the sys- 
tem. In fact, this is another path for the determination 
of q that in the case of unimodal maps provides the same 



results [7] as the alternative methods mentioned above. 
Moreover, rigorous analytical results have recently been 
found for such systems along this line [8]. The tempo- 
ral evolution of the Sq family, allowing to determine qc, 
has also been investigated for diverse other dynamical 
systems [9-12]. The outcome qc 1 has usually been 
interpreted as a signal of nonextensivity. 

Under the scope of the generalized entropies Sq^ we will 
study in the present paper the temporal evolution of a 
cyclic Lotka-Volterra model in the lattice (LLV) [13,14]. 
This model is particularly rich as it exhibits interesting 
features, such as stationary states with spatial patterns 
[13] and fractality [15], which makes it a potential candi- 
date for the applicability of the entropies Sq. The LLV 
is one of the descendants of the original versions con- 
structed by Lotka [16] and Volterra [17] to model auto- 
catalytic chemical reactions and the prey-predator dy- 
namics, respectively, and further adapted to many other 
situations from active transport by proteins [18] to social 
processes [19]. The generalized LLV extends the origi- 
nal scheme to A/" species Ai {i = 1,...,A/') such that, 
cyclically, Ai is "predator" of Ai_i and simultaneously 
"prey" of ^i+i (with ^aT+i = ^i)- Furthermore, the 
dynamics takes place over a lattice, an ingredient that 
introduces new interesting spatial features in comparison 
to the spatially homogeneous mean- field description [20] . 

The criterion of finite entropy rate has been applied be- 
fore to the LLV model in one and two dimensions [12] and 
the relation Qc = 1 — 1/D was conjectured for arbitrary 
D. Now we will extend that study to higher dimensions 
in order to test the conjecture. Moreover, we will essay 
an interpretation of the connection between qc and space 
dimensionality, showing that 7^ 1 is not necessarily 
related to the fractal properties of the dynamics. The 
remaining of the paper is organized as follows. In Sect. 
II we describe in detail the specific LLV model consid- 
ered. In Sect. Ill the time behavior of the entropies Sq 
associated to the LLV is studied. Sect. IV exhibits the 
behavior of Sq under simple transformations of a prob- 
ability density leading to reanalyze in Sect. V the LLV 
dynamics. Final remarks are presented in Sect. VI. 



1 



II. THE GENERALIZED LOTKA-VOLTERRA 
MODEL IN A LATTICE 

Particles of the M different species Ai are localized at 
the sites of a D-dimensional hypercubic lattice. Reac- 
tions between particles of different species occur in bi- 
molecular autocatalytic steps following the scheme 

^2A,+i, (2) 

for i = 1, . . . ,A/', being Atv'+i = ^i, and where < /c^ < 
1 are the kinetic rates. No sites are empty but one of the 
species could be interpreted as representing empty sites 
in the lattice [12,14,15,21]. The dynamics is implemented 
by means of a Monte Carlo (MC) algorithm following the 
details in Refs. [12,14,15]. Basically, at every microscopic 
step: (i) one lattice site is randomly chosen; (ii) one of 
its nearest neighbors is randomly chosen; (iii) if the first 
site is Ai and the neighbor A^+i, the first site changes to 
Aij^i with probability /c^, in accord with scheme (2), oth- 
erwise the site remains unchanged. Each MC step (MCS, 
our unit of time) consists m. N = microscopic steps 
defined above, where L is the linear size of the lattice, so 
that at each MCS all sites are revisited once in average. 
Boundary conditions are periodic. In this paper we will 
deal with three species only {M = 3) and we will focus on 
the symmetric case where all the kinetic rates are equal 
to one (/ci = /c2 = /ca = 1). Moreover, we will restrict our 
study to a particular mode of the dynamics, the "nucleus 
growth" one [12,21]. The lattice is set with all the sites 
filled with species A^ and a "droplet" or nucleus repre- 
sented by a sublattice of small linear size A (i.e., A <C 
is introduced. The droplet contains equal amounts of 
particles of all the three species randomly and uniformly 
distributed. 

As the system evolves according to the MC dynamics, 
the droplet grows and acquires spontaneously a peculiar 
structure [22]. Typical snapshots of the dynamics are 
exhibited in Fig. 1 for lattices in D = 2, 3 and 4 dimen- 
sions. Rings of alternating species develop, repeating the 
sequence Ai^A2^As towards the center of the droplet. 
This pattern is possible since all the kinetic constants 
are equal, then layers have almost the same radial ve- 
locity. The thickness of the rings decreases towards the 
center and thickness fluctuations destroy the most im- 
mersed rings. This behavior has already been observed 
before by Provata and Tsekouras for the two dimensional 
case [21]. They also observed that the destroyed rings 
give rise to a spatial organization of the species in do- 
mains with fractal boundaries typical of the steady state 
in fully occupied periodical lattices [12,15,21]. As the 
dimensionality increases, spatial features remain qualita- 
tively similar but length scales become shorter. In the 
mean- field limit D ^ oo, homogeneity is expected. We 
are going to inspect immediately the evolution of spatial 
patterns from the viewpoint of the generalized entropies 



t = 40 MCS t = 60 MCS t = 80 MCS 




FIG. 1. Snapshots of the dynamics at different times indi- 
cated on the top of the figure, for = 2, 3 and 4. The initial 
condition is a droplet, where all the species are present in 
equal amounts and homogeneously distributed, over a back- 
ground of A3. In all cases the linear size of the lattice is 
L = 100 and the initial nucleus has A = 6. For D > 2 the 
snapshots correspond to sections parallel to one of the hy- 
percube faces and passing by the position of the center of 
the initial droplet. Species are Ai (gray), A2 (black) and A3 
(white) . 



III. TEMPORAL EVOLUTION OF Sq IN THE LLV 

Following previous work [12], we study the temporal 
evolution of the entropies associated to one of the species, 
e.g., Ai. This particular choice does not substantially af- 
fect the results. Nonover lapping windows or sublattices 
{Wi, 1 < i < M} of edge length / covering all the lattice 
are considered. Thus the number of windows must be 
M = (L/l)^ ^ where / is a divisor of L. We associate to 
each window i a probability Pi{t) of being occupied by 
species Ai at time t by counting the number of particles 
of that species ni (i , t) . Then Pi{t) = ni{i^t)/ni {t) , where 
ni{t) is the total number of particles Ai on the lattice. 
The resulting set of probabilities is used to calculate the 
entropy of the lattice, given by the discrete version of 
Eq. (1) (where we have set k = 1) 



The LLV in D = 1 and 2 dimensions has already 
been studied from this viewpoint before [12]. Now we 
will investigate higher dimensional lattices. Fig. 2 shows 
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the generalized entropies as a function of time for di- 
verse values of when the lattice has D = 3 and 
4 dimensions. The results shown in Fig. 2 are not 
qualitatively affected by changing / and L, as long as 
1 < I L. They mainly affect the saturation level, 
gsat. ^ - 1]/(1 - q) for a uniform distribu- 

tion. 
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FIG. 2. Time evolution of the generalized entropies. Sq vs. 
t for various values of q and L — 148 in D = 3 (a), L = 60 
in D = 4 (b). Sq/L vs. t/L for different lattice sizes L with: 
q 2/3 when D = 3 (c), q 3/4 when D = 4 (d). Symbols 
correspond to a single representative numerical experiment. 
In all cases, the window of the partitioning is / = 4 and the 
initial droplet size A = 4. 

Notice in Figs. 2. a and 2.b that as q increases the 
concavity passes from positive to negative, before satu- 
ration, which occurs when the droplet radius becomes of 
the order of the linear size of the lattice. Entropies with 
q yielding constant slope (null concavity) are represented 
as a function of time in Figs. 2.c and 2.d. In all cases, 
constant slope occurs at a value of the entropic index qc 
that as function of D follows the law 



(4) 



in agreement with the relation numerically found for 
D = 1 and 2, and conjectured for generic D previously 
[12]. See Fig. 3. 




FIG. 3. Characteristic value qc as a function of the lattice 
dimensionality D. The dotted straight line corresponds to 
qc = 1 — 1/D. Symbols are the result of numerical experi- 
ments. Errors are of the order of symbol size. For D — 3 and 
4: qc was numerically determined from the data in Fig. 2 as 
the value of q yielding a linear slope. For D = 1 and 2, qc was 
extracted from Ref. [12]. 



IV. BEHAVIOR OF Sq UNDER SIMPLE 
TRANSFORMATIONS OF A PROBABILITY 
DENSITY 



The characteristic entropic index qc is related to the 
number of dimensions D though a simple law. This leads 
to think that a simple mechanism could be behind. Since 
we are dealing with a growth process, we will investigate 
in this Section the relation between qc and D resulting 
from some basic growth mechanisms and connect it with 
Eq. (4). 

Let us analyze the behavior of the entropies Sq under a 
rescaling transformation of an arbitrary probability den- 
sity being x a point in a /^-dimensional space [23]. 
The entropy of the rescaled function pa{x) = p{x / a) / ^ 
where a is a linear length, becomes 



(5) 



As <j ^ oo, the distribution broadens and the entropy 
grows due to the loss of order. If the scaling or stretching 
parameter increases exponentially with time, linear en- 
tropy increase occurs for g = 1. If the scaling parameter 
follows the law a ~ t^/^ with 7 > 0, then the generalized 
entropies increase with time as follows: For q < 1, they 
scale with time as Sq{t) - t'^^-Q)Di/2^ - In t in the 

marginal case ^ = 1, and for q > 1, saturation occurs at 
Sq{oo) = l/{q — 1). Hence a linear increase is achieved 
in the regime g < 1 for 



1 



2 



(6) 



As a relevant particular case let us mention the normal 
diffusive spreading of a D-dimensional Gaussian distribu- 
tion where the stretching parameter is cr^(t) = 2Qt, with 
Q the positive diffusion constant. In this case 7 = 1 then 
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qc = 1 — 2/D. The general case given by Eq. (6) may 
be due to anomalous diffusive motion. A still simpler 
example of rescaling transformation that will be useful 
later consists in a density being non-null and uniform in- 
side a D-dimensional hypersphere and zero outside, such 
that the radius (j{t) grows as a ~ t^/^. Notice that the 
expression for qc of the growing LLV droplet, given by 
Eq. (4), is obtained when the temporal variation of the 
scaling length a is ballistic (7 = 2). 

Now, let us consider another simple process. M ^ N 
windows cover a lattice, like for the LLV in the previous 
Section, with probabilities 



(pi,P2,...,PAr,0,0, ...,0), 



M-N 



l<i<N 



Pi 



1. 



Imagine that windows with non-null probability are repli- 
cated by an integer factor m such that the new set of 
probabilities is 



m ' 



Pi P2 



P2 

m 



Pn 
m 



Pn 
m 



,0,0,...,0) 



M- 



If the replication factor is m ~ (where a is a typical 
linear size of the D-dimensional system) and additionally 
if cr ~ t^/^, then the entropy will be produced in such a 
way that a growth linear in time occurs for qc given by Eq. 
(6). In fact, if a continuous distribution is discretized, the 
associated rescaling process can be thought as a particu- 
lar case of the replication one here considered, a special 
case where replicas are spatially ordered. Nevertheless 
for discrete probabilities, the entropic form Sq does not 
depend on the spatial localization of the windows. 

Summarizing, while for an ordinary exponential 
growth process, the standard entropy Si increases lin- 
early with t; for a process where growth follows a power- 
law in time, Sq has the property of finite asymptotic en- 
tropy rate for some q = qc I. This value of the entropic 
index can result from simple transformations and there- 
fore may be trivially related to the lattice dimensionality. 
These ideas lead us to review the results of the precedent 
Section to see whether simple mechanisms are also lurk- 
ing there. 



interface, we measure A^^, the number of reactive sites 
constituting the interface, and Nf^ the total number of 
reactive faces in inter facial sites. Reactive faces are those 
separating two nearest neighboring sites occupied by dif- 
ferent species. A suitable quantity for our model is 



Nf 



(7) 



where c is the connectivity, i.e, the number of first nearest 
neighbors per site (c = 2D in hypercubic lattices). The 
quantity r represents an averaged measure of the degree 
of reactivity of an inter facial site, and must be also con- 
nected to the roughness of the interface. The reactivity 
r is plotted in Fig. 4 . a as a function of time. It soon 
reaches a stationary value within small fluctuations. For 
comparison, the figure also exhibits the value of r for a 
regularly filled hypersphere with the same total number 
of cells as in the LLV nucleus at each given t. The higher 
the lattice dimensionality, the larger the relative differ- 
ence between the values of r for the two models. The 
roughness of the propagation front in the ID dynamics 
has been studied in detail by Provata and Tsekouras [21]. 
Also in this case, it was shown that after a brief transient 
the rough profile remains stable in average. 
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V. SOME DETAILS OF THE LLV DYNAMICS 

Let us inspect first the roughness of the interface be- 
cause the reaction rate and therefore the dynamics of the 
propagating front are related to that quantity. Since we 
are interested in the interface, one can perform numer- 
ical simulations considering just two species. Hence we 
follow the evolution of an initial nucleus containing only 
Ai particles dropped over a background of A3. Its prop- 
agation is like in the LLV case where a forefront of Ai 
particles governs the spreading of the droplet. Since the 
reaction rate must be proportional to the "extent" of the 



FIG. 4. Propagation of a droplet of Ai in a lattice filled 
with A3, for different values of D indicated in the figure, (a) 
Reactivity r as a function of time for the two-species LLV (full 
lines); for comparison, r for a regularly filled hypersphere with 
the same total number of cells as the corresponding LLV nu- 
cleus is also plotted (dotted lines), (b) Total number of sites 
ni occupied by Ai vs. time. The initial nucleus has linear 
size A = 4. 

If the reactivity remains constant in time, the front 
moves at constant radial speed: R{t) = R{0) -\-vt^ where 

R is an effective radius R cx and the velocity 

V = v{r) an increasing function of its argument for a 
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given D. Although the fractal properties of the inter- 
face are embodied in the velocity, if the roughness soon 
reaches its steady value, then fractality does not affect 
the temporal law of the front propagation. In effect, the 
linear dependence of R with time is observed in numerical 
experiments (see Fig. 4.b). Therefore, following the con- 
siderations of the precedent Section, Eq. (6) with 7 = 2 
must clearly hold in this binary case. 

In the three-species case, although the reaction scheme 
is cyclic, symmetry is broken by the initial condition. 
Species plays the special role of a background what 
in turn makes Ai play the special role of the forefront 
species. After a short transient and before the limits 
of the lattice are reached by the nucleus, the stationar- 
ity of the interface reactivity r, leads to a linear growth 
of the effective radius with time, as in the two species 
case. Consequently, the total number of cells in the nu- 
cleus, protected by the most external ring constituted by 
Ai cells, will increase as . Behind the first ring, the 
three species are equivalent, then the total number of 
cells occupied by each species in the nucleus is expected 
to increase with the same law too. In fact, this be- 
havior is observed in numerical simulations, as shown 
in Fig. 5, although species Ai yields a slightly larger 
slope than A2. Linear growth implies a regime where 
the concentrations in the nucleus are conserved. When 
full occupation of the lattice is attained, the concentra- 
tion of each species q = n^/L^ {i = 1, . . . , 3) fluctuates 
around the stable center predicted by the mean-field the- 
ory: Ci = ki/ Y,j = 1/3, for ah i [14]. 




FIG. 5. Time evolution of the number of cells m and n2 
occupied by species Ai (a) and A2 (b), respectively, for the 
values of the lattice dimensionality indicated in the figure and 
for various lattice sizes L: 150 and 200 (D = 2), 100 and 150 
{D = 3), and 60 and 100 {D = 4), represented by dotted 
and full lines respectively. Time is measured in MCS. Curves 
correspond to single runs. The initial droplet size is A = 6 for 
D = 2 and A 4 for i:> = 3, 4. 



We have seen that the production of the species in the 
nucleus is such that after a short transient their con- 
centrations remain approximately constant. Thus, in 
the extreme case when windows have the size of a cell 
(/ = 1), a linear increase of Sq with time occurs clearly 
for q given by Eq. (4). In fact there will be ni windows 
with non-null probability 1/ni, where ni{t) ~ t^, then 
Sq{t) = (ni"^-l)/(l-^), so that: Sq - t, if {l-q)D = 1. 
In the opposite extreme where windows are so large that 
all non-empty windows have approximately the same oc- 
cupation number, that same temporal dependence is ob- 
tained too. 
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FIG. 6. Cumulative distribution of occupation numbers at 
different times (in MCS) indicated in the figure for D = 2 (a) 
and D — (b). is the number of windows occupied by 

less than i cells of species A\ (1 < i < l^). In all cases, the 
window edges have length / = 4. ¥01 D ^ 2, (L, A) = (200, 6), 
and for D = 3, (L, A) = (148, 4). Insets: multiplication factor 
Ctt/Ctt^a^ as a function of time (in MCS), where the horizon- 
tal bars mean average over all the occupation numbers at a 
given t, and tmax = 160 (a), 120 (b). 

Let us analyze what happens for intermediate window 
sizes. Clusters of cells of each species like in fully occu- 
pied lattices [12,15] appear in the interior of the nucleus. 
It is not counterintuitive the idea that as the volume of 
the nucleus is increasing, the size distribution of the ag- 
glomerates becomes stationary after a transient. To test 
this idea one can calculate for instance the occupation 
numbers of the windows, intimately related to entropy 
computation, and see how their distribution evolves in 
time. Fig. 6 exhibits the number of windows ^{i) occu- 
pied by at most i cells of species Ai (1 < i < /^), for 
different time instants. The occupation number for z = 
is dismissed whereas it does not contribute to the entropy 
and a cumulative probability is considered in order to get 
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smoother curves. Notice that the curves in the log-hnear 
plot are practicahy parahel, indicating that they differ in 
a multiphcative factor. Moreover, the rephcation factor 
increases fohowing the law as illustrated in the inset 
of Fig. 6 and in accord with previous considerations. All 
this means that a replication of the windows as described 
in Sect. IV with 7 = 2 is going on. 

In conclusion, although the underlying dynamics is 
complex, once the properties of spatial structures attain 
steady values, growth is controlled by simple laws. More 
specifically, the resulting growth processes are power-law 
in time with scaling exponents trivially related to the lat- 
tice dimensionality (with cr ^ t or, equivalently, 7 = 2). 
Then, the linear increase of Sq with time for qc given by 
Eq. (4) is expected. 



VI. FINAL REMARKS 

In this paper we employed the criterion of finite asymp- 
totic growth rate of Sq to determine Qc in the D- 
dimensional LLV. We verified for D = 3 and D = 4 a 
conjecture for the connection between Qc and D previ- 
ously proposed [12], namely, Qc = 1 — 1/ D. Moreover we 
exhibited the mechanisms leading to such relation. 

The box-counting method applied to the 2D LLV [15] 
has shown that the boundaries of the small domains of a 
given species are approximately fractal [15] with a frac- 
tal dimension df that depends on the reaction rates ki. 
Analogous results are expected in higher dimensions. We 
have seen that, for the present model and for the partic- 
ular probability assignment considered, Sq entropy pro- 
duction does not capture directly the fractality of spatial 
patterns. Of course, since the fractal dimension df must 
depend on the lattice spatial dimension the character- 
istic entropic index Qc (a function of D) results in some 
way connected to df. But qc is not determined by the 
degree of fractality, it is only determined by indepen- 
dently of the nature (fractal or not) of the growing core. 
This is so because the properties of spatial patterns, such 
as the interface roughness, soon reach steady values. 

In a process where the number of occupied cells in- 
creases exponentially with time the usual BGS entropy 
Si increases linearly with t. For a growth process that 
is not exponential, one can not expect a linear increase 
of the standard entropy Si. If growth occurs following 
a power-law in time, Sq has the property of finite pro- 
duction rate for some q = Qc 1. Then 7^ 1 is the 
expected behavior in theses cases and it is not necessar- 
ily connected to the complex features of a system. So qc 
may result trivially related to the lattice dimensionality. 

We have shown that the stationarity of certain prop- 
erties of the LLV dynamics determines a growth process 
linear in time (7 = 2), yielding relation (4). Although 
simple, this relation is in some way a consequence of the 
complex LLV dynamics. As perspectives, one can not 
exclude the possibility that in other dynamical regimes 



of the LLV, the increase of Sq can reflect complex fea- 
tures directly. Also, it could be insightful to review 
previous works in the literature by taking into consid- 
eration the present results. This might be especially 
fruitful in cases where characteristic indexes of the form 
qc = 1 — 2/{jD) have also been found, as in the interest- 
ing study of Galilean- invariant lattice Boltzmann models 
of fluids [11]. 
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